home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / HFIELD.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-05-02  |  45.7 KB  |  2,278 lines
  1. /****************************************************************************
  2. *                   hfield.c
  3. *
  4. *    This file implements the height field shape primitive.  The shape is
  5. *    implemented as a collection of triangles which are calculated as
  6. *    needed.  The basic intersection routine first computes the rays
  7. *    intersection with the box marking the limits of the shape, then
  8. *    follows the line from one intersection point to the other, testing the
  9. *    two triangles which form the pixel for an intersection with the ray at
  10. *    each step.
  11. *
  12. *    height field added by Doug Muir with lots of advice and support
  13. *    from David Buck and Drew Wells.
  14. *
  15. *  from Persistence of Vision(tm) Ray Tracer
  16. *  Copyright 1996 Persistence of Vision Team
  17. *---------------------------------------------------------------------------
  18. *  NOTICE: This source code file is provided so that users may experiment
  19. *  with enhancements to POV-Ray and to port the software to platforms other
  20. *  than those supported by the POV-Ray Team.  There are strict rules under
  21. *  which you are permitted to use this file.  The rules are in the file
  22. *  named POVLEGAL.DOC which should be distributed with this file. If
  23. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  24. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  25. *  Forum.  The latest version of POV-Ray may be found there as well.
  26. *
  27. * This program is based on the popular DKB raytracer version 2.12.
  28. * DKBTrace was originally written by David K. Buck.
  29. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  30. *
  31. *****************************************************************************/
  32.  
  33. /****************************************************************************
  34. *
  35. *  Explanation:
  36. *
  37. *    -
  38. *
  39. *  Syntax:
  40. *
  41. *  ---
  42. *
  43. *  Aug 1994 : Merged functions for CSG height fields into functions for
  44. *             non-CSG height fiels. Moved all height field map related
  45. *             data into one data structure. Fixed memory problems. [DB]
  46. *
  47. *  Feb 1995 : Major rewrite of the height field intersection tests. [DB]
  48. *
  49. *****************************************************************************/
  50.  
  51. #include "frame.h"
  52. #include "povray.h"
  53. #include "vector.h"
  54. #include "povproto.h"
  55. #include "hfield.h"
  56. #include "matrices.h"
  57. #include "objects.h"
  58.  
  59.  
  60.  
  61. /*****************************************************************************
  62. * Local preprocessor defines
  63. ******************************************************************************/
  64.  
  65. #define sign(x) (((x) >= 0.0) ? 1.0 : -1.0)
  66.  
  67. #define Get_Height(x, z, HField) ((DBL)(HField)->Data->Map[(z)][(x)])
  68.  
  69. /* Small offest. */
  70.  
  71. #define HFIELD_OFFSET 0.001
  72.  
  73.  
  74.  
  75. /*****************************************************************************
  76. * Static functions
  77. ******************************************************************************/
  78.  
  79. static DBL normalize PARAMS((VECTOR A, VECTOR B));
  80.  
  81. static DBL stretch PARAMS((DBL x));
  82.  
  83. static void smooth_height_field PARAMS((HFIELD *HField, int xsize, int zsize));
  84.  
  85. static int intersect_pixel PARAMS((int x,int z,RAY *Ray,HFIELD *HField,DBL height1,DBL height2));
  86.  
  87. static int add_single_normal PARAMS((HF_VAL **data, int xsize, int zsize,
  88.   int x0, int z0,int x1, int z1,int x2, int z2,VECTOR N));
  89.  
  90. static int  All_HField_Intersections PARAMS((OBJECT *Object,RAY *Ray,ISTACK *Depth_Stack));
  91. static int  Inside_HField PARAMS((VECTOR IPoint,OBJECT *Object));
  92. static void HField_Normal PARAMS((VECTOR Result,OBJECT *Object,INTERSECTION *Inter));
  93. static void Translate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  94. static void Rotate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  95. static void Scale_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  96. static void Invert_HField PARAMS((OBJECT *Object));
  97. static void Transform_HField PARAMS((OBJECT *Object,TRANSFORM *Trans));
  98. static void *Copy_HField PARAMS((OBJECT *Object));
  99. static void Destroy_HField PARAMS((OBJECT *Object));
  100.  
  101. static int dda_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start, HFIELD_BLOCK *Block));
  102. static int block_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start));
  103. static void build_hfield_blocks PARAMS((HFIELD *HField));
  104.  
  105.  
  106.  
  107. /*****************************************************************************
  108. * Local variables
  109. ******************************************************************************/
  110.  
  111. METHODS HField_Methods =
  112. {
  113.   All_HField_Intersections,
  114.   Inside_HField, HField_Normal,
  115.   Copy_HField, Translate_HField, Rotate_HField,
  116.   Scale_HField, Transform_HField, Invert_HField, Destroy_HField
  117. };
  118.  
  119. static ISTACK *HField_Stack;
  120. static RAY *RRay;
  121. static DBL mindist, maxdist;
  122.  
  123.  
  124. /*****************************************************************************
  125. *
  126. * FUNCTION
  127. *
  128. *   All_HField_Intersections
  129. *
  130. * INPUT
  131. *   
  132. * OUTPUT
  133. *   
  134. * RETURNS
  135. *   
  136. * AUTHOR
  137. *
  138. *   Doug Muir, David Buck, Drew Wells
  139. *   
  140. * DESCRIPTION
  141. *
  142. *   -
  143. *
  144. * CHANGES
  145. *
  146. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  147. *
  148. ******************************************************************************/
  149.  
  150. static int All_HField_Intersections(Object, Ray, Depth_Stack)
  151. OBJECT *Object;
  152. RAY *Ray;
  153. ISTACK *Depth_Stack;
  154. {
  155.   int Side1, Side2;
  156.   VECTOR Start;
  157.   RAY Temp_Ray;
  158.   DBL depth1, depth2;
  159.   HFIELD *HField = (HFIELD *) Object;
  160.  
  161.   Increase_Counter(stats[Ray_HField_Tests]);
  162.  
  163.   MInvTransPoint(Temp_Ray.Initial, Ray->Initial, HField->Trans);
  164.   MInvTransDirection(Temp_Ray.Direction, Ray->Direction, HField->Trans);
  165.  
  166. #ifdef HFIELD_EXTRA_STATS
  167.   Increase_Counter(stats[Ray_HField_Box_Tests]);
  168. #endif
  169.  
  170.   if (!Intersect_Box(&Temp_Ray,HField->bounding_box,&depth1,&depth2,&Side1,&Side2))
  171.   {
  172.     return(FALSE);
  173.   }
  174.  
  175. #ifdef HFIELD_EXTRA_STATS
  176.   Increase_Counter(stats[Ray_HField_Box_Tests_Succeeded]);
  177. #endif
  178.  
  179.   HField->Data->cache_pos = 0;
  180.  
  181.   if (depth1 < EPSILON)
  182.   {
  183.     depth1 = EPSILON;
  184.   }
  185.  
  186.   VEvaluateRay(Start, Temp_Ray.Initial, depth1, Temp_Ray.Direction);
  187.  
  188.   mindist = depth1;
  189.   maxdist = depth2;
  190.  
  191.   HField_Stack = Depth_Stack;
  192.  
  193.   RRay = Ray;
  194.  
  195.   if (block_traversal(&Temp_Ray, HField, Start))
  196.   {
  197.     Increase_Counter(stats[Ray_HField_Tests_Succeeded]);
  198.  
  199.     return(TRUE);
  200.   }
  201.   else
  202.   {
  203.     return(FALSE);
  204.   }
  205. }
  206.  
  207.  
  208.  
  209. /*****************************************************************************
  210. *
  211. * FUNCTION
  212. *
  213. *   Inside_HField
  214. *
  215. * INPUT
  216. *   
  217. * OUTPUT
  218. *   
  219. * RETURNS
  220. *   
  221. * AUTHOR
  222. *
  223. *   Doug Muir, David Buck, Drew Wells
  224. *   
  225. * DESCRIPTION
  226. *
  227. *   -
  228. *
  229. * CHANGES
  230. *
  231. *   -
  232. *
  233. ******************************************************************************/
  234.  
  235. static int Inside_HField (IPoint, Object)
  236. VECTOR IPoint;
  237. OBJECT *Object;
  238. {
  239.   HFIELD *HField = (HFIELD *) Object;
  240.   int px, pz;
  241.   DBL x,z,y1,y2,y3,water, dot1, dot2;
  242.   VECTOR Local_Origin, H_Normal, Test;
  243.  
  244.   MInvTransPoint(Test, IPoint, HField->Trans);
  245.  
  246.   water = HField->bounding_box->bounds[0][Y];
  247.  
  248.   if ((Test[X] < 0.0) || (Test[X] >= HField->bounding_box->bounds[1][X]) ||
  249.       (Test[Z] < 0.0) || (Test[Z] >= HField->bounding_box->bounds[1][Z]))
  250.   {
  251.     return(Test_Flag(HField, INVERTED_FLAG));
  252.   }
  253.  
  254.   if (Test[Y] >= HField->bounding_box->bounds[1][Y])
  255.   {
  256.     return(Test_Flag(HField, INVERTED_FLAG));
  257.   }
  258.  
  259.   if (Test[Y] < water)
  260.   {
  261.     return(!Test_Flag(HField, INVERTED_FLAG));
  262.   }
  263.  
  264.   px = (int)Test[X];
  265.   pz = (int)Test[Z];
  266.  
  267.   x = Test[X] - (DBL)px;
  268.   z = Test[Z] - (DBL)pz;
  269.  
  270.   if ((x+z) < 1.0)
  271.   {
  272.     y1 = max(Get_Height(px,   pz,   HField), water);
  273.     y2 = max(Get_Height(px+1, pz,   HField), water);
  274.     y3 = max(Get_Height(px,   pz+1, HField), water);
  275.  
  276.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  277.  
  278.     Make_Vector(H_Normal, y1-y2, 1.0, y1-y3);
  279.   }
  280.   else
  281.   {
  282.     px = (int)ceil(Test[X]);
  283.     pz = (int)ceil(Test[Z]);
  284.  
  285.     y1 = max(Get_Height(px,   pz,   HField), water);
  286.     y2 = max(Get_Height(px-1, pz,   HField), water);
  287.     y3 = max(Get_Height(px,   pz-1, HField), water);
  288.  
  289.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  290.  
  291.     Make_Vector(H_Normal, y2-y1, 1.0, y3-y1);
  292.   }
  293.  
  294.   VDot(dot1, Test, H_Normal);
  295.   VDot(dot2, Local_Origin, H_Normal);
  296.  
  297.   if (dot1 < dot2)
  298.   {
  299.     return(!Test_Flag(HField, INVERTED_FLAG));
  300.   }
  301.  
  302.   return(Test_Flag(HField, INVERTED_FLAG));
  303. }
  304.  
  305.  
  306.  
  307. /*****************************************************************************
  308. *
  309. * FUNCTION
  310. *
  311. *   HField_Normal
  312. *
  313. * INPUT
  314. *   
  315. * OUTPUT
  316. *   
  317. * RETURNS
  318. *   
  319. * AUTHOR
  320. *
  321. *   Doug Muir, David Buck, Drew Wells
  322. *   
  323. * DESCRIPTION
  324. *
  325. *   -
  326. *
  327. * CHANGES
  328. *
  329. *   -
  330. *
  331. ******************************************************************************/
  332.  
  333. static void HField_Normal (Result, Object, Inter)
  334. OBJECT *Object;
  335. VECTOR Result;
  336. INTERSECTION *Inter;
  337. {
  338.   HFIELD *HField = (HFIELD *) Object;
  339.   int px,pz, i;
  340.   DBL x,z,y1,y2,y3,u,v;
  341.   VECTOR Local_Origin;
  342.   VECTOR n[5];
  343.  
  344.   MInvTransPoint(Local_Origin, Inter->IPoint, HField->Trans);
  345.  
  346.   for (i = 0; i < HField->Data->cache_pos; i++)
  347.   {
  348.     if ((Local_Origin[X] == HField->Data->Normal_Cache[i].fx) &&
  349.        (Local_Origin[Z] == HField->Data->Normal_Cache[i].fz))
  350.     {
  351.       Assign_Vector(Result,HField->Data->Normal_Cache[i].normal);
  352.  
  353.       MTransNormal(Result,Result,HField->Trans);
  354.  
  355.       VNormalize(Result,Result);
  356.  
  357.       return;
  358.     }
  359.   }
  360.  
  361.   px = (int)Local_Origin[X];
  362.   pz = (int)Local_Origin[Z];
  363.  
  364.   x = Local_Origin[X] - (DBL)px;
  365.   z = Local_Origin[Z] - (DBL)pz;
  366.  
  367.   if (Test_Flag(HField, SMOOTHED_FLAG))
  368.   {
  369.     n[0][X] = HField->Data->Normals[pz][px][0];
  370.     n[0][Y] = HField->Data->Normals[pz][px][1];
  371.     n[0][Z] = HField->Data->Normals[pz][px][2];
  372.     n[1][X] = HField->Data->Normals[pz][px+1][0];
  373.     n[1][Y] = HField->Data->Normals[pz][px+1][1];
  374.     n[1][Z] = HField->Data->Normals[pz][px+1][2];
  375.     n[2][X] = HField->Data->Normals[pz+1][px][0];
  376.     n[2][Y] = HField->Data->Normals[pz+1][px][1];
  377.     n[2][Z] = HField->Data->Normals[pz+1][px][2];
  378.     n[3][X] = HField->Data->Normals[pz+1][px+1][0];
  379.     n[3][Y] = HField->Data->Normals[pz+1][px+1][1];
  380.     n[3][Z] = HField->Data->Normals[pz+1][px+1][2];
  381.  
  382.     x = stretch(x);
  383.     z = stretch(z);
  384.  
  385.     u = (1.0 - x);
  386.     v = (1.0 - z);
  387.  
  388.     Result[X] = v*(u*n[0][X] + x*n[1][X]) + z*(u*n[2][X] + x*n[3][X]);
  389.     Result[Y] = v*(u*n[0][Y] + x*n[1][Y]) + z*(u*n[2][Y] + x*n[3][Y]);
  390.     Result[Z] = v*(u*n[0][Z] + x*n[1][Z]) + z*(u*n[2][Z] + x*n[3][Z]);
  391.   }
  392.   else
  393.   {
  394.     if ((x+z) <= 1.0)
  395.     {
  396.       /* Lower triangle. */
  397.  
  398.       y1 = Get_Height(px,   pz,   HField);
  399.       y2 = Get_Height(px+1, pz,   HField);
  400.       y3 = Get_Height(px,   pz+1, HField);
  401.  
  402.       Make_Vector(Result, y1-y2, 1.0, y1-y3);
  403.     }
  404.     else
  405.     {
  406.       /* Upper triangle. */
  407.  
  408.       y1 = Get_Height(px+1, pz+1, HField);
  409.       y2 = Get_Height(px,   pz+1, HField);
  410.       y3 = Get_Height(px+1, pz,   HField);
  411.  
  412.       Make_Vector(Result, y2-y1, 1.0, y3-y1);
  413.     }
  414.   }
  415.  
  416.   MTransNormal(Result, Result, HField->Trans);
  417.  
  418.   VNormalize(Result, Result);
  419. }
  420.  
  421.  
  422.  
  423. /*****************************************************************************
  424. *
  425. * FUNCTION
  426. *
  427. *   stretch
  428. *
  429. * INPUT
  430. *   
  431. * OUTPUT
  432. *   
  433. * RETURNS
  434. *   
  435. * AUTHOR
  436. *
  437. *   Doug Muir, David Buck, Drew Wells
  438. *   
  439. * DESCRIPTION
  440. *
  441. *   -
  442. *
  443. * CHANGES
  444. *
  445. *   -
  446. *
  447. ******************************************************************************/
  448.  
  449. static DBL stretch (x)
  450. DBL x;
  451. {
  452.   if (x <= 0.5)
  453.   {
  454.     x = 2.0 * x * x;
  455.   }
  456.   else
  457.   {
  458.     x = 1.0 - (2.0 * (1.0 - x) * (1.0 - x));
  459.   }
  460.  
  461.   return(x);
  462. }
  463.  
  464.  
  465.  
  466. /*****************************************************************************
  467. *
  468. * FUNCTION
  469. *
  470. *   normalize
  471. *
  472. * INPUT
  473. *   
  474. * OUTPUT
  475. *   
  476. * RETURNS
  477. *   
  478. * AUTHOR
  479. *
  480. *   Doug Muir, David Buck, Drew Wells
  481. *   
  482. * DESCRIPTION
  483. *
  484. *   -
  485. *
  486. * CHANGES
  487. *
  488. *   -
  489. *
  490. ******************************************************************************/
  491.  
  492. static DBL normalize(A, B)
  493. VECTOR A, B;
  494. {
  495.   VLength(VTemp, B);
  496.  
  497.   if (fabs(VTemp) > EPSILON)
  498.   {
  499.     VInverseScale(A, B, VTemp);
  500.   }
  501.   else
  502.   {
  503.     Make_Vector(A, 0.0, 1.0, 0.0);
  504.   }
  505.  
  506.   return(VTemp);
  507. }
  508.  
  509.  
  510.  
  511. /*****************************************************************************
  512. *
  513. * FUNCTION
  514. *
  515. *   intersect_pixel
  516. *
  517. * INPUT
  518. *   
  519. * OUTPUT
  520. *   
  521. * RETURNS
  522. *   
  523. * AUTHOR
  524. *
  525. *   Doug Muir, David Buck, Drew Wells
  526. *   
  527. * DESCRIPTION
  528. *
  529. *   -
  530. *
  531. * CHANGES
  532. *
  533. *   -
  534. *
  535. ******************************************************************************/
  536.  
  537. static int intersect_pixel(x, z, Ray, HField, height1, height2)
  538. int x;
  539. int z;
  540. RAY *Ray;
  541. HFIELD *HField;
  542. DBL height1, height2;
  543. {
  544.   int Found;
  545.   DBL dot, depth1, depth2;
  546.   DBL s, t, y1, y2, y3, y4;
  547.   DBL min_y2_y3, max_y2_y3;
  548.   DBL max_height, min_height;
  549.   VECTOR P, N, V1;
  550.  
  551. #ifdef HFIELD_EXTRA_STATS
  552.   Increase_Counter(stats[Ray_HField_Cell_Tests]);
  553. #endif
  554.  
  555.   y1 = Get_Height(x,   z,   HField);
  556.   y2 = Get_Height(x+1, z,   HField);
  557.   y3 = Get_Height(x,   z+1, HField);
  558.   y4 = Get_Height(x+1, z+1, HField);
  559.  
  560.   /* Do we hit this cell at all? */
  561.  
  562.   if (y2 < y3)
  563.   {
  564.     min_y2_y3 = y2;
  565.     max_y2_y3 = y3;
  566.   }
  567.   else
  568.   {
  569.     min_y2_y3 = y3;
  570.     max_y2_y3 = y2;
  571.   }
  572.  
  573.   min_height = min(min(y1, y4), min_y2_y3);
  574.   max_height = max(max(y1, y4), max_y2_y3);
  575.  
  576.   if ((max_height < height1) || (min_height > height2))
  577.   {
  578.     return(FALSE);
  579.   }
  580.  
  581. #ifdef HFIELD_EXTRA_STATS
  582.   Increase_Counter(stats[Ray_HField_Cell_Tests_Succeeded]);
  583. #endif
  584.  
  585.   Found = FALSE;
  586.  
  587.   /* Check if we'll hit first triangle. */
  588.  
  589.   min_height = min(y1, min_y2_y3);
  590.   max_height = max(y1, max_y2_y3);
  591.  
  592.   if ((max_height >= height1) && (min_height <= height2))
  593.   {
  594. #ifdef HFIELD_EXTRA_STATS
  595.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  596. #endif
  597.  
  598.     /* Set up triangle. */
  599.  
  600.     Make_Vector(P, (DBL)x, y1, (DBL)z);
  601.  
  602.     /*
  603.      * Calculate the normal vector from:
  604.      *
  605.      * N = V2 x V1, with V1 = <1, y2-y1, 0>, V2 = <0, y3-y1, 1>.
  606.      */
  607.  
  608.     Make_Vector(N, y1-y2, 1.0, y1-y3);
  609.  
  610.     /* Now intersect the triangle. */
  611.  
  612.     VDot(dot, N, Ray->Direction);
  613.  
  614.     if ((dot > EPSILON) || (dot < -EPSILON))
  615.     {
  616.       VSub(V1, P, Ray->Initial);
  617.  
  618.       VDot(depth1, N, V1);
  619.  
  620.       depth1 /= dot;
  621.  
  622.       if ((depth1 >= mindist) && (depth1 <= maxdist))
  623.       {
  624.         s = Ray->Initial[X] + depth1 * Ray->Direction[X] - (DBL)x;
  625.         t = Ray->Initial[Z] + depth1 * Ray->Direction[Z] - (DBL)z;
  626.  
  627.         if ((s >= -0.0001) && (t >= -0.0001) && ((s+t) <= 1.0001))
  628.         {
  629. #ifdef HFIELD_EXTRA_STATS
  630.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  631. #endif
  632.  
  633.           VEvaluateRay(P, RRay->Initial, depth1, RRay->Direction);
  634.  
  635.           if (Point_In_Clip(P, HField->Clip))
  636.           {
  637.             push_entry(depth1, P, (OBJECT *)HField, HField_Stack);
  638.  
  639.             Found = TRUE;
  640.  
  641.             /* Cache normal. */
  642.  
  643.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  644.             {
  645.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  646.               {
  647.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  648.  
  649.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  650.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  651.  
  652.                 HField->Data->cache_pos++;
  653.               }
  654.             }
  655.           }
  656.         }
  657.       }
  658.     }
  659.   }
  660.  
  661.   /* Check if we'll hit second triangle. */
  662.  
  663.   min_height = min(y4, min_y2_y3);
  664.   max_height = max(y4, max_y2_y3);
  665.  
  666.   if ((max_height >= height1) && (min_height <= height2))
  667.   {
  668. #ifdef HFIELD_EXTRA_STATS
  669.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  670. #endif
  671.  
  672.     /* Set up triangle. */
  673.  
  674.     Make_Vector(P, (DBL)(x+1), y4, (DBL)(z+1));
  675.  
  676.     /*
  677.      * Calculate the normal vector from:
  678.      *
  679.      * N = V2 x V1, with V1 = <-1, y3-y4, 0>, V2 = <0, y2-y4, -1>.
  680.      */
  681.  
  682.     Make_Vector(N, y3-y4, 1.0, y2-y4);
  683.  
  684.     /* Now intersect the triangle. */
  685.  
  686.     VDot(dot, N, Ray->Direction);
  687.  
  688.     if ((dot > EPSILON) || (dot < -EPSILON))
  689.     {
  690.       VSub(V1, P, Ray->Initial);
  691.  
  692.       VDot(depth2, N, V1);
  693.  
  694.       depth2 /= dot;
  695.  
  696.       if ((depth2 >= mindist) && (depth2 <= maxdist))
  697.       {
  698.         s = Ray->Initial[X] + depth2 * Ray->Direction[X] - (DBL)x;
  699.         t = Ray->Initial[Z] + depth2 * Ray->Direction[Z] - (DBL)z;
  700.  
  701.         if ((s <= 1.0001) && (t <= 1.0001) && ((s+t) >= 0.9999))
  702.         {
  703. #ifdef HFIELD_EXTRA_STATS
  704.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  705. #endif
  706.  
  707.           VEvaluateRay(P, RRay->Initial, depth2, RRay->Direction);
  708.  
  709.           if (Point_In_Clip(P, HField->Clip))
  710.           {
  711.             push_entry(depth2, P, (OBJECT *)HField, HField_Stack);
  712.  
  713.             Found = TRUE;
  714.  
  715.             /* Cache normal. */
  716.  
  717.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  718.             {
  719.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  720.               {
  721.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  722.  
  723.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  724.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  725.  
  726.                 HField->Data->cache_pos++;
  727.               }
  728.             }
  729.           }
  730.         }
  731.       }
  732.     }
  733.   }
  734.  
  735.   return(Found);
  736. }
  737.  
  738.  
  739.  
  740. /*****************************************************************************
  741. *
  742. * FUNCTION
  743. *
  744. *   add_single_normal
  745. *
  746. * INPUT
  747. *   
  748. * OUTPUT
  749. *   
  750. * RETURNS
  751. *   
  752. * AUTHOR
  753. *
  754. *   Doug Muir, David Buck, Drew Wells
  755. *   
  756. * DESCRIPTION
  757. *
  758. *   -
  759. *
  760. * CHANGES
  761. *
  762. *   -
  763. *
  764. ******************************************************************************/
  765.  
  766. static int add_single_normal(data, xsize, zsize, x0, z0, x1, z1, x2, z2, N)
  767. HF_VAL **data;
  768. int xsize,zsize,x0,z0,x1,z1,x2,z2;
  769. VECTOR N;
  770. {
  771.   VECTOR v0, v1, v2, t0, t1, Nt;
  772.  
  773.   if ((x0 < 0) || (z0 < 0) ||
  774.       (x1 < 0) || (z1 < 0) ||
  775.       (x2 < 0) || (z2 < 0) ||
  776.       (x0 > xsize) || (z0 > zsize) ||
  777.       (x1 > xsize) || (z1 > zsize) ||
  778.       (x2 > xsize) || (z2 > zsize))
  779.   {
  780.     return(0);
  781.   }
  782.   else
  783.   {
  784.     Make_Vector(v0, x0, (DBL)data[z0][x0], z0);
  785.     Make_Vector(v1, x1, (DBL)data[z1][x1], z1);
  786.     Make_Vector(v2, x2, (DBL)data[z2][x2], z2);
  787.  
  788.     VSub(t0, v2, v0);
  789.     VSub(t1, v1, v0);
  790.  
  791.     VCross(Nt, t0, t1);
  792.  
  793.     normalize(Nt, Nt);
  794.  
  795.     if (Nt[Y] < 0.0)
  796.     {
  797.       VScaleEq(Nt, -1.0);
  798.     }
  799.  
  800.     VAddEq(N, Nt);
  801.  
  802.     return(1);
  803.   }
  804. }
  805.  
  806.  
  807.  
  808. /*****************************************************************************
  809. *
  810. * FUNCTION
  811. *
  812. *   smooth_height_field
  813. *
  814. * INPUT
  815. *   
  816. * OUTPUT
  817. *   
  818. * RETURNS
  819. *   
  820. * AUTHOR
  821. *
  822. *   Doug Muir, David Buck, Drew Wells
  823. *   
  824. * DESCRIPTION
  825. *
  826. *   Given a height field that only contains an elevation grid, this
  827. *   routine will walk through the data and produce averaged normals
  828. *   for all points on the grid.
  829. *
  830. * CHANGES
  831. *
  832. *   -
  833. *
  834. ******************************************************************************/
  835.  
  836. static void smooth_height_field(HField, xsize, zsize)
  837. HFIELD *HField;
  838. int xsize;
  839. int zsize;
  840. {
  841.   int i, j, k;
  842.   VECTOR N;
  843.   HF_VAL **map = HField->Data->Map;
  844.  
  845.   /* First off, allocate all the memory needed to store the normal information */
  846.  
  847.   HField->Data->Normals_Height = zsize+1;
  848.  
  849.   HField->Data->Normals = (HF_Normals **)POV_MALLOC((zsize+1)*sizeof(HF_Normals *), "height field normals");
  850.  
  851.   for (i = 0; i <= zsize; i++)
  852.   {
  853.     HField->Data->Normals[i] = (HF_Normals *)POV_MALLOC((xsize+1)*sizeof(HF_Normals), "height field normals");
  854.   }
  855.  
  856.   /*
  857.    * For now we will do it the hard way - by generating the normals
  858.    * individually for each elevation point.
  859.    */
  860.  
  861.   for (i = 0; i < zsize; i++)
  862.   {
  863.     COOPERATE_1
  864.  
  865.     for (j = 0; j < xsize; j++)
  866.     {
  867.       Make_Vector(N, 0.0, 0.0, 0.0);
  868.  
  869.       k = 0;
  870.  
  871.       k += add_single_normal(map, xsize, zsize, j, i, j+1, i, j, i+1, N);
  872.       k += add_single_normal(map, xsize, zsize, j, i, j, i+1, j-1, i, N);
  873.       k += add_single_normal(map, xsize, zsize, j, i, j-1, i, j, i-1, N);
  874.       k += add_single_normal(map, xsize, zsize, j, i, j, i-1, j+1, i, N);
  875.  
  876.       if (k == 0)
  877.       {
  878.         Error ("Failed to find any normals at: (%d, %d).\n", i, j);
  879.       }
  880.  
  881.       normalize(N, N);
  882.  
  883.       HField->Data->Normals[i][j][0] = (short)(32767 * N[X]);
  884.       HField->Data->Normals[i][j][1] = (short)(32767 * N[Y]);
  885.       HField->Data->Normals[i][j][2] = (short)(32767 * N[Z]);
  886.     }
  887.   }
  888. }
  889.  
  890.  
  891.  
  892. /*****************************************************************************
  893. *
  894. * FUNCTION
  895. *
  896. *   Compute_HField
  897. *
  898. * INPUT
  899. *
  900. * OUTPUT
  901. *
  902. * RETURNS
  903. *
  904. * AUTHOR
  905. *
  906. *   Doug Muir, David Buck, Drew Wells
  907. *
  908. * DESCRIPTION
  909. *
  910. *   Copy image data into height field map. Create bounding blocks
  911. *   for the block traversal. Calculate normals for smoothed height fields.
  912. *
  913. * CHANGES
  914. *
  915. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  916. *
  917. ******************************************************************************/
  918.  
  919. void Compute_HField(HField, Image)
  920. HFIELD *HField;
  921. IMAGE *Image;
  922. {
  923.   int x, z, max_x, max_z;
  924.   int temp1 = 0, temp2 = 0;
  925.   HF_VAL min_y, max_y, temp_y;
  926.  
  927.   /* Get height field map size. */
  928.  
  929.   max_x = Image->iwidth;
  930.  
  931.   if (Image->File_Type == POT_FILE)
  932.   {
  933.     max_x = max_x / 2;
  934.   }
  935.  
  936.   max_z = Image->iheight;
  937.  
  938.   /* Allocate memory for map. */
  939.  
  940.   HField->Data->Map = (HF_VAL **)POV_MALLOC(max_z*sizeof(HF_VAL *), "height field");
  941.  
  942.   for (z = 0; z < max_z; z++)
  943.   {
  944.     HField->Data->Map[z] = (HF_VAL *)POV_MALLOC(max_x*sizeof(HF_VAL), "height field");
  945.   }
  946.  
  947.   /* Copy map. */
  948.  
  949.   min_y = 65535L;
  950.   max_y = 0;
  951.  
  952.   for (z = 0; z < max_z; z++)
  953.   {
  954.     COOPERATE_1
  955.     for (x = 0; x < max_x; x++)
  956.     {
  957.  
  958.       switch (Image->File_Type)
  959.       {
  960.         case GIF_FILE:
  961.  
  962.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  963.           temp2 = 0;
  964.  
  965.           break;
  966.  
  967.         case POT_FILE:
  968.  
  969.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  970.           temp2 = Image->data.map_lines[max_z - z - 1][x + max_x];
  971.  
  972.           break;
  973.  
  974.  
  975.         case PPM_FILE:
  976.  
  977.           temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  978.           temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  979.  
  980.           break;
  981.  
  982.         case PGM_FILE:
  983.         case TGA_FILE:
  984.         case PNG_FILE:
  985.  
  986.           if (Image->Colour_Map == NULL)
  987.           {
  988.             temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  989.             temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  990.           }
  991.           else
  992.           {
  993.             temp1 = Image->data.map_lines[max_z - z - 1][x];
  994.             temp2 = 0;
  995.           }
  996.  
  997.           break;
  998.  
  999.         default:
  1000.  
  1001.           Error("Unknown image type in Compute_HField().\n");
  1002.       }
  1003.  
  1004.       temp_y = (HF_VAL)(256*temp1 + temp2);
  1005.  
  1006.       HField->Data->Map[z][x] = temp_y;
  1007.  
  1008.       min_y = min(min_y, temp_y);
  1009.       max_y = max(max_y, temp_y);
  1010.     }
  1011.   }
  1012.  
  1013.   /* Resize bounding box. */
  1014.  
  1015.   HField->Data->min_y = min_y;
  1016.   HField->Data->max_y = max_y;
  1017.  
  1018.   HField->bounding_box->bounds[0][Y] = max((DBL)min_y, HField->bounding_box->bounds[0][Y]) - HFIELD_OFFSET;
  1019.   HField->bounding_box->bounds[1][Y] = (DBL)max_y + HFIELD_OFFSET;
  1020.  
  1021.   /* Compute smoothed height field. */
  1022.  
  1023.   if (Test_Flag(HField, SMOOTHED_FLAG))
  1024.   {
  1025.     smooth_height_field(HField, max_x-1, max_z-1);
  1026.   }
  1027.  
  1028.   HField->Data->max_x = max_x-2;
  1029.   HField->Data->max_z = max_z-2;
  1030.  
  1031.   build_hfield_blocks(HField);
  1032. }
  1033.  
  1034.  
  1035.  
  1036. /*****************************************************************************
  1037. *
  1038. * FUNCTION
  1039. *
  1040. *   build_hfield_blocks
  1041. *
  1042. * INPUT
  1043. *   
  1044. * OUTPUT
  1045. *   
  1046. * RETURNS
  1047. *   
  1048. * AUTHOR
  1049. *
  1050. *   Dieter Bayer
  1051. *   
  1052. * DESCRIPTION
  1053. *
  1054. *   Create the bounding block hierarchy used by the block traversal.
  1055. *
  1056. * CHANGES
  1057. *
  1058. *   Feb 1995 : Creation.
  1059. *
  1060. ******************************************************************************/
  1061.  
  1062. static void build_hfield_blocks(HField)
  1063. HFIELD *HField;
  1064. {
  1065.   int x, z, nx, nz, wx, wz;
  1066.   int i, j;
  1067.   int xmin, xmax, zmin, zmax;
  1068.   DBL y, ymin, ymax, water;
  1069.  
  1070.   /* Get block size. */
  1071.  
  1072.   nx = max(1, (int)(sqrt((DBL)HField->Data->max_x)));
  1073.   nz = max(1, (int)(sqrt((DBL)HField->Data->max_z)));
  1074.  
  1075.   /* Get dimensions of sub-block. */
  1076.  
  1077.   wx = (int)ceil((DBL)(HField->Data->max_x + 2) / (DBL)nx);
  1078.   wz = (int)ceil((DBL)(HField->Data->max_z + 2) / (DBL)nz);
  1079.  
  1080.   /* Increase number of sub-blocks if necessary. */
  1081.  
  1082.   if (nx * wx < HField->Data->max_x + 2)
  1083.   {
  1084.     nx++;
  1085.   }
  1086.  
  1087.   if (nz * wz < HField->Data->max_z + 2)
  1088.   {
  1089.     nz++;
  1090.   }
  1091.  
  1092.   if (!Test_Flag(HField, HIERARCHY_FLAG) || ((nx == 1) && (nz == 1)))
  1093.   {
  1094.     /* We don't want a bounding hierarchy. Just use one block. */
  1095.  
  1096.     HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(sizeof(HFIELD_BLOCK *), "height field blocks");
  1097.  
  1098.     HField->Data->Block[0] = (HFIELD_BLOCK *)POV_MALLOC(sizeof(HFIELD_BLOCK), "height field blocks");
  1099.  
  1100.     HField->Data->Block[0][0].xmin = 0;
  1101.     HField->Data->Block[0][0].xmax = HField->Data->max_x;
  1102.     HField->Data->Block[0][0].zmin = 0;
  1103.     HField->Data->Block[0][0].zmax = HField->Data->max_z;
  1104.  
  1105.     HField->Data->Block[0][0].ymin = HField->bounding_box->bounds[0][Y];
  1106.     HField->Data->Block[0][0].ymax = HField->bounding_box->bounds[1][Y];
  1107.  
  1108.     HField->Data->block_max_x = 1;
  1109.     HField->Data->block_max_z = 1;
  1110.  
  1111.     HField->Data->block_width_x = HField->Data->max_x + 2;
  1112.     HField->Data->block_width_z = HField->Data->max_y + 2;
  1113.  
  1114. /*
  1115.     Debug_Info("\nHeight field: %d x %d (1 x 1 blocks)", HField->Data->max_x+2, HField->Data->max_z+2);
  1116. */
  1117.  
  1118.     return;
  1119.   }
  1120.  
  1121. /*
  1122.   Debug_Info("\nHeight field: %d x %d (%d x %d blocks)", HField->Data->max_x+2, HField->Data->max_z+2, nx, nz);
  1123. */
  1124.  
  1125.   /* Allocate memory for blocks. */
  1126.  
  1127.   HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(nz*sizeof(HFIELD_BLOCK *), "height field blocks");
  1128.  
  1129.   /* Store block information. */
  1130.  
  1131.   HField->Data->block_max_x = nx;
  1132.   HField->Data->block_max_z = nz;
  1133.  
  1134.   HField->Data->block_width_x = wx;
  1135.   HField->Data->block_width_z = wz;
  1136.  
  1137.   water = HField->bounding_box->bounds[0][Y];
  1138.  
  1139.   for (z = 0; z < nz; z++)
  1140.   {
  1141.     COOPERATE_1
  1142.  
  1143.     HField->Data->Block[z] = (HFIELD_BLOCK *)POV_MALLOC(nx*sizeof(HFIELD_BLOCK), "height field blocks");
  1144.  
  1145.     for (x = 0; x < nx; x++)
  1146.     {
  1147.       /* Get block's borders. */
  1148.  
  1149.       xmin = x * wx;
  1150.       zmin = z * wz;
  1151.  
  1152.       xmax = min((x + 1) * wx - 1, HField->Data->max_x);
  1153.       zmax = min((z + 1) * wz - 1, HField->Data->max_z);
  1154.  
  1155.       /* Find min. and max. height in current block. */
  1156.  
  1157.       ymin = BOUND_HUGE;
  1158.       ymax = -BOUND_HUGE;
  1159.  
  1160.       for (i = xmin; i <= xmax+1; i++)
  1161.       {
  1162.         for (j = zmin; j <= zmax+1; j++)
  1163.         {
  1164.           y = Get_Height(i, j, HField);
  1165.  
  1166.           ymin = min(ymin, y);
  1167.           ymax = max(ymax, y);
  1168.         }
  1169.       }
  1170.  
  1171.       /* Store block's borders. */
  1172.  
  1173.       HField->Data->Block[z][x].xmin = xmin;
  1174.       HField->Data->Block[z][x].xmax = xmax;
  1175.       HField->Data->Block[z][x].zmin = zmin;
  1176.       HField->Data->Block[z][x].zmax = zmax;
  1177.  
  1178.       HField->Data->Block[z][x].ymin = max(ymin, water) - HFIELD_OFFSET;
  1179.       HField->Data->Block[z][x].ymax = ymax + HFIELD_OFFSET;
  1180.     }
  1181.   }
  1182. }
  1183.  
  1184.  
  1185.  
  1186. /*****************************************************************************
  1187. *
  1188. * FUNCTION
  1189. *
  1190. *   Translate_HField
  1191. *
  1192. * INPUT
  1193. *   
  1194. * OUTPUT
  1195. *   
  1196. * RETURNS
  1197. *   
  1198. * AUTHOR
  1199. *
  1200. *   Doug Muir, David Buck, Drew Wells
  1201. *
  1202. * DESCRIPTION
  1203. *
  1204. *   -
  1205. *
  1206. * CHANGES
  1207. *
  1208. *   -
  1209. *
  1210. ******************************************************************************/
  1211.  
  1212. static void Translate_HField (Object, Vector, Trans)
  1213. OBJECT *Object;
  1214. VECTOR Vector;
  1215. TRANSFORM *Trans;
  1216. {
  1217.   Transform_HField(Object, Trans);
  1218. }
  1219.  
  1220.  
  1221.  
  1222. /*****************************************************************************
  1223. *
  1224. * FUNCTION
  1225. *
  1226. *   Rotate_HField
  1227. *
  1228. * INPUT
  1229. *
  1230. * OUTPUT
  1231. *
  1232. * RETURNS
  1233. *
  1234. * AUTHOR
  1235. *
  1236. *   Doug Muir, David Buck, Drew Wells
  1237. *
  1238. * DESCRIPTION
  1239. *
  1240. *   -
  1241. *
  1242. * CHANGES
  1243. *
  1244. *   -
  1245. *
  1246. ******************************************************************************/
  1247.  
  1248. static void Rotate_HField (Object, Vector, Trans)
  1249. OBJECT *Object;
  1250. VECTOR Vector;
  1251. TRANSFORM *Trans;
  1252. {
  1253.   Transform_HField(Object, Trans);
  1254. }
  1255.  
  1256.  
  1257.  
  1258. /*****************************************************************************
  1259. *
  1260. * FUNCTION
  1261. *
  1262. *   Scale_HField
  1263. *
  1264. * INPUT
  1265. *   
  1266. * OUTPUT
  1267. *   
  1268. * RETURNS
  1269. *   
  1270. * AUTHOR
  1271. *
  1272. *   Doug Muir, David Buck, Drew Wells
  1273. *   
  1274. * DESCRIPTION
  1275. *
  1276. *   -
  1277. *
  1278. * CHANGES
  1279. *
  1280. *   -
  1281. *
  1282. ******************************************************************************/
  1283.  
  1284. static void Scale_HField (Object, Vector, Trans)
  1285. OBJECT *Object;
  1286. VECTOR Vector;
  1287. TRANSFORM *Trans;
  1288. {
  1289.   Transform_HField(Object, Trans);
  1290. }
  1291.  
  1292.  
  1293.  
  1294. /*****************************************************************************
  1295. *
  1296. * FUNCTION
  1297. *
  1298. *   Invert_HField
  1299. *
  1300. * INPUT
  1301. *   
  1302. * OUTPUT
  1303. *   
  1304. * RETURNS
  1305. *   
  1306. * AUTHOR
  1307. *
  1308. *   Doug Muir, David Buck, Drew Wells
  1309. *   
  1310. * DESCRIPTION
  1311. *
  1312. *   -
  1313. *
  1314. * CHANGES
  1315. *
  1316. *   -
  1317. *
  1318. ******************************************************************************/
  1319.  
  1320. static void Invert_HField (Object)
  1321. OBJECT *Object;
  1322. {
  1323.   Invert_Flag(Object, INVERTED_FLAG);
  1324. }
  1325.  
  1326.  
  1327.  
  1328. /*****************************************************************************
  1329. *
  1330. * FUNCTION
  1331. *
  1332. *   Transform_HField
  1333. *
  1334. * INPUT
  1335. *   
  1336. * OUTPUT
  1337. *   
  1338. * RETURNS
  1339. *   
  1340. * AUTHOR
  1341. *
  1342. *   Doug Muir, David Buck, Drew Wells
  1343. *   
  1344. * DESCRIPTION
  1345. *
  1346. *   -
  1347. *
  1348. * CHANGES
  1349. *
  1350. *   -
  1351. *
  1352. ******************************************************************************/
  1353.  
  1354. static void Transform_HField (Object, Trans)
  1355. OBJECT *Object;
  1356. TRANSFORM *Trans;
  1357. {
  1358.   Compose_Transforms(((HFIELD *)Object)->Trans, Trans);
  1359.  
  1360.   Compute_HField_BBox((HFIELD *)Object);
  1361. }
  1362.  
  1363.  
  1364.  
  1365. /*****************************************************************************
  1366. *
  1367. * FUNCTION
  1368. *
  1369. *   Create_HField
  1370. *
  1371. * INPUT
  1372. *
  1373. * OUTPUT
  1374. *   
  1375. * RETURNS
  1376. *   
  1377. * AUTHOR
  1378. *
  1379. *   Doug Muir, David Buck, Drew Wells
  1380. *   
  1381. * DESCRIPTION
  1382. *
  1383. *   Allocate and intialize a height field.
  1384. *
  1385. * CHANGES
  1386. *
  1387. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1388. *
  1389. ******************************************************************************/
  1390.  
  1391. HFIELD *Create_HField()
  1392. {
  1393.   HFIELD *New;
  1394.  
  1395.   /* Allocate height field. */
  1396.  
  1397.   New = (HFIELD *)POV_MALLOC(sizeof(HFIELD), "height field");
  1398.  
  1399.   INIT_OBJECT_FIELDS(New, HFIELD_OBJECT, &HField_Methods)
  1400.  
  1401.   /* Always uses Trans so always create one. */
  1402.  
  1403.   New->Trans = Create_Transform();
  1404.  
  1405.   New->bounding_box = Create_Box();
  1406.  
  1407.   /* Allocate height field data. */
  1408.  
  1409.   New->Data = (HFIELD_DATA *)POV_MALLOC(sizeof(HFIELD_DATA), "height field");
  1410.  
  1411.   New->Data->References = 1;
  1412.  
  1413.   New->Data->cache_pos = 0;
  1414.  
  1415.   New->Data->Normals_Height = 0;
  1416.  
  1417.   New->Data->Map     = NULL;
  1418.   New->Data->Normals = NULL;
  1419.  
  1420.   New->Data->max_x = 0;
  1421.   New->Data->max_z = 0;
  1422.  
  1423.   New->Data->block_max_x = 0;
  1424.   New->Data->block_max_z = 0;
  1425.  
  1426.   New->Data->block_width_x = 0;
  1427.   New->Data->block_width_z = 0;
  1428.  
  1429.   Set_Flag(New, HIERARCHY_FLAG);
  1430.  
  1431.   return(New);
  1432. }
  1433.  
  1434.  
  1435.  
  1436. /*****************************************************************************
  1437. *
  1438. * FUNCTION
  1439. *
  1440. *   Copy_HField
  1441. *
  1442. * INPUT
  1443. *   
  1444. * OUTPUT
  1445. *   
  1446. * RETURNS
  1447. *   
  1448. * AUTHOR
  1449. *
  1450. *   Doug Muir, David Buck, Drew Wells
  1451. *   
  1452. * DESCRIPTION
  1453. *
  1454. *   NOTE: The height field data is not copied, only the number of references
  1455. *         is counted, so that Destray_HField() knows if it can be destroyed.
  1456. *
  1457. * CHANGES
  1458. *
  1459. *   -
  1460. *
  1461. ******************************************************************************/
  1462.  
  1463. static void *Copy_HField(Object)
  1464. OBJECT *Object;
  1465. {
  1466.   HFIELD *New;
  1467.  
  1468.   New = Create_HField();
  1469.  
  1470.   /* Destroy obsolete things created in Create_HField(). */
  1471.  
  1472.   Destroy_Transform(New->Trans);
  1473.  
  1474.   Destroy_Box((OBJECT *)(New->bounding_box));
  1475.  
  1476.   POV_FREE (New->Data);
  1477.  
  1478.   /* Copy height field. */
  1479.  
  1480.   *New = *((HFIELD *)Object);
  1481.  
  1482.   New->Trans = Copy_Transform(((HFIELD *)Object)->Trans);
  1483.  
  1484.   New->bounding_box = Copy_Box((OBJECT *)(((HFIELD *)Object)->bounding_box));
  1485.  
  1486.   New->Data->References++;
  1487.  
  1488.   return(New);
  1489. }
  1490.  
  1491.  
  1492.  
  1493. /*****************************************************************************
  1494. *
  1495. * FUNCTION
  1496. *
  1497. *   Destroy_HField
  1498. *
  1499. * INPUT
  1500. *   
  1501. * OUTPUT
  1502. *   
  1503. * RETURNS
  1504. *   
  1505. * AUTHOR
  1506. *
  1507. *   Doug Muir, David Buck, Drew Wells
  1508. *   
  1509. * DESCRIPTION
  1510. *
  1511. *   NOTE: The height field data is destroyed if it's no longer
  1512. *         used by any copy.
  1513. *
  1514. * CHANGES
  1515. *
  1516. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1517. *
  1518. ******************************************************************************/
  1519.  
  1520. static void Destroy_HField (Object)
  1521. OBJECT *Object;
  1522. {
  1523.   int i;
  1524.   HFIELD *HField = (HFIELD *)Object;
  1525.  
  1526.   Destroy_Transform(HField->Trans);
  1527.  
  1528.   Destroy_Box((OBJECT *)(HField->bounding_box));
  1529.  
  1530.   if (--(HField->Data->References) == 0)
  1531.   {
  1532.     if (HField->Data->Map != NULL)
  1533.     {
  1534.       for (i = 0; i < HField->Data->max_z+2; i++)
  1535.       {
  1536.         if (HField->Data->Map[i] != NULL)
  1537.         {
  1538.           POV_FREE (HField->Data->Map[i]);
  1539.         }
  1540.       }
  1541.  
  1542.       POV_FREE (HField->Data->Map);
  1543.     }
  1544.  
  1545.     if (HField->Data->Normals != NULL)
  1546.     {
  1547.       for (i = 0; i < HField->Data->Normals_Height; i++)
  1548.       {
  1549.         POV_FREE (HField->Data->Normals[i]);
  1550.       }
  1551.  
  1552.       POV_FREE (HField->Data->Normals);
  1553.     }
  1554.  
  1555.     if (HField->Data->Block != NULL)
  1556.     {
  1557.       for (i = 0; i < HField->Data->block_max_z; i++)
  1558.       {
  1559.         POV_FREE(HField->Data->Block[i]);
  1560.       }
  1561.  
  1562.       POV_FREE(HField->Data->Block);
  1563.     }
  1564.  
  1565.     POV_FREE (HField->Data);
  1566.   }
  1567.  
  1568.   POV_FREE (Object);
  1569. }
  1570.  
  1571.  
  1572.  
  1573. /*****************************************************************************
  1574. *
  1575. * FUNCTION
  1576. *
  1577. *   Compute_HField_BBox
  1578. *
  1579. * INPUT
  1580. *
  1581. *   HField - Height field
  1582. *   
  1583. * OUTPUT
  1584. *
  1585. *   HField
  1586. *   
  1587. * RETURNS
  1588. *   
  1589. * AUTHOR
  1590. *
  1591. *   Dieter Bayer
  1592. *   
  1593. * DESCRIPTION
  1594. *
  1595. *   Calculate the bounding box of a height field.
  1596. *
  1597. * CHANGES
  1598. *
  1599. *   Aug 1994 : Creation.
  1600. *
  1601. ******************************************************************************/
  1602.  
  1603. void Compute_HField_BBox(HField)
  1604. HFIELD *HField;
  1605. {
  1606.   Assign_BBox_Vect(HField->BBox.Lower_Left, HField->bounding_box->bounds[0]);
  1607.  
  1608.   VSub (HField->BBox.Lengths, HField->bounding_box->bounds[1], HField->bounding_box->bounds[0]);
  1609.  
  1610.   if (HField->Trans != NULL)
  1611.   {
  1612.     Recompute_BBox(&HField->BBox, HField->Trans);
  1613.   }
  1614. }
  1615.  
  1616.  
  1617.  
  1618. /*****************************************************************************
  1619. *
  1620. * FUNCTION
  1621. *
  1622. *   dda_traversal
  1623. *
  1624. * INPUT
  1625. *
  1626. *   Ray    - Current ray
  1627. *   HField - Height field
  1628. *   Start  - Start point for the walk
  1629. *   Block  - Sub-block of the height field to traverse
  1630. *   
  1631. * OUTPUT
  1632. *   
  1633. * RETURNS
  1634. *
  1635. *   int - TRUE if intersection was found
  1636. *   
  1637. * AUTHOR
  1638. *
  1639. *   Dieter Bayer
  1640. *   
  1641. * DESCRIPTION
  1642. *
  1643. *   Traverse the grid cell of the height field using a modified DDA.
  1644. *
  1645. *   Based on the following article:
  1646. *
  1647. *     Musgrave, F. Kenton, "Grid Tracing: Fast Ray Tracing for Height
  1648. *     Fields", Research Report YALEU-DCS-RR-39, Yale University, July 1988
  1649. *
  1650. *   You should note that there are (n-1) x (m-1) grid cells in a height
  1651. *   field of (image) size n x m. A grid cell (i,j), 0 <= i <= n-1,
  1652. *   0 <= j <= m-1, extends from x = i to x = i + 1 - epsilon and
  1653. *   y = j to y = j + 1 -epsilon.
  1654. *
  1655. * CHANGES
  1656. *
  1657. *   Feb 1995 : Creation.
  1658. *
  1659. ******************************************************************************/
  1660.  
  1661. static int dda_traversal(Ray, HField, Start, Block)
  1662. RAY *Ray;
  1663. HFIELD *HField;
  1664. VECTOR Start;
  1665. HFIELD_BLOCK *Block;
  1666. {
  1667.   int found;
  1668.   int xmin, xmax, zmin, zmax;
  1669.   int x, z, signx, signz;
  1670.   DBL ymin, ymax, y1, y2;
  1671.   DBL px, pz, dx, dy, dz;
  1672.   DBL delta, error, x0, z0;
  1673.   DBL neary, fary, deltay;
  1674.  
  1675.   /* Setup DDA. */
  1676.  
  1677.   found = FALSE;
  1678.  
  1679.   px = Start[X];
  1680.   pz = Start[Z];
  1681.  
  1682.   /* Get dimensions of current block. */
  1683.  
  1684.   xmin = Block->xmin;
  1685.   xmax = min(Block->xmax + 1, HField->Data->max_x);
  1686.   zmin = Block->zmin;
  1687.   zmax = min(Block->zmax + 1, HField->Data->max_z);
  1688.  
  1689.   ymin = min(Start[Y], Block->ymin) - EPSILON;
  1690.   ymax = max(Start[Y], Block->ymax) + EPSILON;
  1691.  
  1692.   /* Check for illegal grid values (caused by numerical inaccuracies). */
  1693.  
  1694.   if (px < (DBL)xmin)
  1695.   {
  1696.     if (px < (DBL)xmin - HFIELD_OFFSET)
  1697.     {
  1698.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1699.  
  1700.       return(FALSE);
  1701.     }
  1702.     else
  1703.     {
  1704.       px = (DBL)xmin;
  1705.     }
  1706.   }
  1707.   else
  1708.   {
  1709.     if (px > (DBL)xmax + 1.0 - EPSILON)
  1710.     {
  1711.       if (px > (DBL)xmax + 1.0 + EPSILON)
  1712.       {
  1713.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1714.  
  1715.         return(FALSE);
  1716.       }
  1717.       else
  1718.       {
  1719.         px = (DBL)xmax + 1.0 - EPSILON;
  1720.       }
  1721.     }
  1722.   }
  1723.  
  1724.   if (pz < (DBL)zmin)
  1725.   {
  1726.     if (pz < (DBL)zmin - HFIELD_OFFSET)
  1727.     {
  1728.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1729.  
  1730.       return(FALSE);
  1731.     }
  1732.     else
  1733.     {
  1734.       pz = (DBL)zmin;
  1735.     }
  1736.   }
  1737.   else
  1738.   {
  1739.     if (pz > (DBL)zmax + 1.0 - EPSILON)
  1740.     {
  1741.       if (pz > (DBL)zmax + 1.0 + EPSILON)
  1742.       {
  1743.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1744.  
  1745.         return(FALSE);
  1746.       }
  1747.       else
  1748.       {
  1749.         pz = (DBL)zmax + 1.0 - EPSILON;
  1750.       }
  1751.     }
  1752.   }
  1753.  
  1754.   dx = Ray->Direction[X];
  1755.   dy = Ray->Direction[Y];
  1756.   dz = Ray->Direction[Z];
  1757.  
  1758.   /*
  1759.    * Here comes the DDA algorithm.
  1760.    */
  1761.  
  1762.   /* Choose algorithm depending on the driving axis. */
  1763.  
  1764.   if (fabs(dx) >= fabs(dz))
  1765.   {
  1766.     /*
  1767.      * X-axis is driving axis.
  1768.      */
  1769.  
  1770.     delta = fabs(dz / dx);
  1771.  
  1772.     x = (int)px;
  1773.     z = (int)pz;
  1774.  
  1775.     x0 = px - floor(px);
  1776.     z0 = pz - floor(pz);
  1777.  
  1778.     signx = sign(dx);
  1779.     signz = sign(dz);
  1780.  
  1781.     /* Get initial error. */
  1782.  
  1783.     if (dx >= 0.0)
  1784.     {
  1785.       if (dz >= 0.0)
  1786.       {
  1787.         error = z0 + delta * (1.0 - x0) - 1.0;
  1788.       }
  1789.       else
  1790.       {
  1791.         error = -(z0 - delta * (1.0 - x0));
  1792.       }
  1793.     }
  1794.     else
  1795.     {
  1796.       if (dz >= 0.0)
  1797.       {
  1798.         error = z0 + delta * x0 - 1.0;
  1799.       }
  1800.       else
  1801.       {
  1802.         error = -(z0 - delta * x0);
  1803.       }
  1804.     }
  1805.  
  1806.     /* Get y differential. */
  1807.  
  1808.     deltay = dy / fabs(dx);
  1809.  
  1810.     if (dx >= 0.0)
  1811.     {
  1812.       neary = Start[Y] - x0 * deltay;
  1813.  
  1814.       fary = neary + deltay;
  1815.     }
  1816.     else
  1817.     {
  1818.       neary = Start[Y] - (1.0 - x0) * deltay;
  1819.  
  1820.       fary = neary + deltay;
  1821.     }
  1822.  
  1823.     /* Step through the cells. */
  1824.  
  1825.     do
  1826.     {
  1827.       if (neary < fary)
  1828.       {
  1829.         y1 = neary;
  1830.         y2 = fary;
  1831.       }
  1832.       else
  1833.       {
  1834.         y1 = fary;
  1835.         y2 = neary;
  1836.       }
  1837.  
  1838.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1839.       {
  1840.         if (HField->Type & IS_CHILD_OBJECT)
  1841.         {
  1842.           found = TRUE;
  1843.         }
  1844.         else
  1845.         {
  1846.           return(TRUE);
  1847.         }
  1848.       }
  1849.  
  1850.       if (error > EPSILON)
  1851.       {
  1852.         z += signz;
  1853.  
  1854.         if ((z < zmin) || (z > zmax))
  1855.         {
  1856.           break;
  1857.         }
  1858.         else
  1859.         {
  1860.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1861.           {
  1862.             if (HField->Type & IS_CHILD_OBJECT)
  1863.             {
  1864.               found = TRUE;
  1865.             }
  1866.             else
  1867.             {
  1868.               return(TRUE);
  1869.             }
  1870.           }
  1871.         }
  1872.  
  1873.         error--;
  1874.       }
  1875.       else
  1876.       {
  1877.         if (error > -EPSILON)
  1878.         {
  1879.           z += signz;
  1880.  
  1881.           error--;
  1882.         }
  1883.       }
  1884.  
  1885.       x += signx;
  1886.  
  1887.       error += delta;
  1888.  
  1889.       neary = fary;
  1890.  
  1891.       fary += deltay;
  1892.     }
  1893.     while ((neary >= ymin) && (neary <= ymax) && (x >= xmin) && (x <= xmax) && (z >= zmin) && (z <= zmax));
  1894.   }
  1895.   else
  1896.   {
  1897.     /*
  1898.      * Z-axis is driving axis.
  1899.      */
  1900.  
  1901.     delta = fabs(dx / dz);
  1902.  
  1903.     x = (int)px;
  1904.     z = (int)pz;
  1905.  
  1906.     x0 = px - floor(px);
  1907.     z0 = pz - floor(pz);
  1908.  
  1909.     signx = sign(dx);
  1910.     signz = sign(dz);
  1911.  
  1912.     /* Get initial error. */
  1913.  
  1914.     if (dz >= 0.0)
  1915.     {
  1916.       if (dx >= 0.0)
  1917.       {
  1918.         error = x0 + delta * (1.0 - z0) - 1.0;
  1919.       }
  1920.       else
  1921.       {
  1922.         error = -(x0 - delta * (1.0 - z0));
  1923.       }
  1924.     }
  1925.     else
  1926.     {
  1927.       if (dx >= 0.0)
  1928.       {
  1929.         error = x0 + delta * z0 - 1.0;
  1930.       }
  1931.       else
  1932.       {
  1933.         error = -(x0 - delta * z0);
  1934.       }
  1935.     }
  1936.  
  1937.     /* Get y differential. */
  1938.  
  1939.     deltay = dy / fabs(dz);
  1940.  
  1941.     if (dz >= 0.0)
  1942.     {
  1943.       neary = Start[Y] - z0 * deltay;
  1944.  
  1945.       fary = neary + deltay;
  1946.     }
  1947.     else
  1948.     {
  1949.       neary = Start[Y] - (1.0 - z0) * deltay;
  1950.  
  1951.       fary = neary + deltay;
  1952.     }
  1953.  
  1954.     /* Step through the cells. */
  1955.  
  1956.     do
  1957.     {
  1958.       if (neary < fary)
  1959.       {
  1960.         y1 = neary;
  1961.         y2 = fary;
  1962.       }
  1963.       else
  1964.       {
  1965.         y1 = fary;
  1966.         y2 = neary;
  1967.       }
  1968.  
  1969.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1970.       {
  1971.         if (HField->Type & IS_CHILD_OBJECT)
  1972.         {
  1973.           found = TRUE;
  1974.         }
  1975.         else
  1976.         {
  1977.           return(TRUE);
  1978.         }
  1979.       }
  1980.  
  1981.       if (error > EPSILON)
  1982.       {
  1983.         x += signx;
  1984.  
  1985.         if ((x < xmin) || (x > xmax))
  1986.         {
  1987.           break;
  1988.         }
  1989.         else
  1990.         {
  1991.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1992.           {
  1993.             if (HField->Type & IS_CHILD_OBJECT)
  1994.             {
  1995.               found = TRUE;
  1996.             }
  1997.             else
  1998.             {
  1999.               return(TRUE);
  2000.             }
  2001.           }
  2002.         }
  2003.  
  2004.         error--;
  2005.       }
  2006.       else
  2007.       {
  2008.         if (error > -EPSILON)
  2009.         {
  2010.           x += signx;
  2011.  
  2012.           error--;
  2013.         }
  2014.       }
  2015.  
  2016.       z += signz;
  2017.  
  2018.       error += delta;
  2019.  
  2020.       neary = fary;
  2021.  
  2022.       fary += deltay;
  2023.     }
  2024.     while ((neary >= ymin-EPSILON) && (neary <= ymax+EPSILON) &&
  2025.            (x >= xmin) && (x <= xmax) &&
  2026.            (z >= zmin) && (z <= zmax));
  2027.   }
  2028.  
  2029.   return(found);
  2030. }
  2031.  
  2032.  
  2033.  
  2034. /*****************************************************************************
  2035. *
  2036. * FUNCTION
  2037. *
  2038. *   block_traversal
  2039. *
  2040. * INPUT
  2041. *
  2042. *   Ray    - Current ray
  2043. *   HField - Height field
  2044. *   Start  - Start point for the walk
  2045. *
  2046. * OUTPUT
  2047. *
  2048. * RETURNS
  2049. *
  2050. *   int - TRUE if intersection was found
  2051. *   
  2052. * AUTHOR
  2053. *
  2054. *   Dieter Bayer
  2055. *
  2056. * DESCRIPTION
  2057. *
  2058. *   Traverse the blocks of the height field.
  2059. *
  2060. * CHANGES
  2061. *
  2062. *   Feb 1995 : Creation.
  2063. *
  2064. ******************************************************************************/
  2065.  
  2066. static int block_traversal(Ray, HField, Start)
  2067. RAY *Ray;
  2068. HFIELD *HField;
  2069. VECTOR Start;
  2070. {
  2071.   int xmax, zmax;
  2072.   int x, z, nx, nz, signx, signz;
  2073.   int found = FALSE;
  2074.   int dx_zero, dz_zero;
  2075.   DBL px, pz, dx, dy, dz;
  2076.   DBL ymin, ymax, y1, y2;
  2077.   DBL neary, fary;
  2078.   DBL k1, k2, dist;
  2079.   VECTOR nearP, farP;
  2080.   HFIELD_BLOCK *Block;
  2081.  
  2082.   px = Start[X];
  2083.   pz = Start[Z];
  2084.  
  2085.   dx = Ray->Direction[X];
  2086.   dy = Ray->Direction[Y];
  2087.   dz = Ray->Direction[Z];
  2088.  
  2089.   /* First test for 'perpendicular' rays. */
  2090.  
  2091.   if ((fabs(dx) < EPSILON) && (fabs(dz) < EPSILON))
  2092.   {
  2093.     x = (int)px;
  2094.     z = (int)pz;
  2095.  
  2096.     neary = Start[Y];
  2097.  
  2098.     if (dy >= 0.0)
  2099.     {
  2100.       fary = 65536.0;
  2101.     }
  2102.     else
  2103.     {
  2104.       fary = 0.0;
  2105.     }
  2106.  
  2107.     return(intersect_pixel(x, z, Ray, HField, min(neary, fary), max(neary, fary)));
  2108.   }
  2109.  
  2110.   /* If we don't have blocks we just step through the grid. */
  2111.  
  2112.   if ((HField->Data->block_max_x <= 1) && (HField->Data->block_max_z <= 1))
  2113.   {
  2114.     return(dda_traversal(Ray, HField, Start, &HField->Data->Block[0][0]));
  2115.   }
  2116.  
  2117.   /* Get dimensions of grid. */
  2118.  
  2119.   xmax = HField->Data->block_max_x;
  2120.   zmax = HField->Data->block_max_z;
  2121.  
  2122.   ymin = (DBL)HField->Data->min_y - EPSILON;
  2123.   ymax = (DBL)HField->Data->max_y + EPSILON;
  2124.  
  2125.   dx_zero = (fabs(dx) < EPSILON);
  2126.   dz_zero = (fabs(dz) < EPSILON);
  2127.  
  2128.   signx = sign(dx);
  2129.   signz = sign(dz);
  2130.  
  2131.   /* Walk on the block grid. */
  2132.  
  2133.   px /= HField->Data->block_width_x;
  2134.   pz /= HField->Data->block_width_z;
  2135.  
  2136.   x = (int)px;
  2137.   z = (int)pz;
  2138.  
  2139.   Assign_Vector(nearP, Start);
  2140.  
  2141.   neary = Start[Y];
  2142.  
  2143.   /*
  2144.    * Here comes the block walk algorithm.
  2145.    */
  2146.  
  2147.   do
  2148.   {
  2149. #ifdef HFIELD_EXTRA_STATS
  2150.     Increase_Counter(stats[Ray_HField_Block_Tests]);
  2151. #endif
  2152.  
  2153.     /* Get current block. */
  2154.  
  2155.     Block = &HField->Data->Block[z][x];
  2156.  
  2157.     /* Intersect ray with bounding planes. */
  2158.  
  2159.     if (dx_zero)
  2160.     {
  2161.       k1 = BOUND_HUGE;
  2162.     }
  2163.     else
  2164.     {
  2165.       if (signx >= 0)
  2166.       {
  2167.         k1 = ((DBL)Block->xmax + 1.0 - Ray->Initial[X]) / dx;
  2168.       }
  2169.       else
  2170.       {
  2171.         k1 = ((DBL)Block->xmin - Ray->Initial[X]) / dx;
  2172.       }
  2173.     }
  2174.  
  2175.     if (dz_zero)
  2176.     {
  2177.       k2 = BOUND_HUGE;
  2178.     }
  2179.     else
  2180.     {
  2181.       if (signz >= 0)
  2182.       {
  2183.         k2 = ((DBL)Block->zmax + 1.0 - Ray->Initial[Z]) / dz;
  2184.       }
  2185.       else
  2186.       {
  2187.         k2 = ((DBL)Block->zmin - Ray->Initial[Z]) / dz;
  2188.       }
  2189.     }
  2190.  
  2191.     /* Figure out the indices of the next block. */
  2192.  
  2193.     if ((k1 < k2 - EPSILON) && (k1 > 0.0))
  2194.     {
  2195.       /* Step along the x-axis. */
  2196.  
  2197.       dist = k1;
  2198.  
  2199.       nx = x + signx;
  2200.       nz = z;
  2201.     }
  2202.     else
  2203.     {
  2204.       if ((k1 < k2 + EPSILON) && (k1 > 0.0))
  2205.       {
  2206.         /* Step along both axis (very rare case). */
  2207.  
  2208.         dist = k1;
  2209.  
  2210.         nx = x + signx;
  2211.         nz = z + signz;
  2212.       }
  2213.       else
  2214.       {
  2215.         /* Step along the z-axis. */
  2216.  
  2217.         dist = k2;
  2218.  
  2219.         nx = x;
  2220.         nz = z + signz;
  2221.       }
  2222.     }
  2223.  
  2224.     /* Get point where ray leaves current block. */
  2225.  
  2226.     VEvaluateRay(farP, Ray->Initial, dist, Ray->Direction);
  2227.  
  2228.     fary = farP[Y];
  2229.  
  2230.     if (neary < fary)
  2231.     {
  2232.       y1 = neary;
  2233.       y2 = fary;
  2234.     }
  2235.     else
  2236.     {
  2237.       y1 = fary;
  2238.       y2 = neary;
  2239.     }
  2240.  
  2241.     /* Can we hit current block at all? */
  2242.  
  2243.     if ((y1 <= (DBL)Block->ymax + EPSILON) && (y2 >= (DBL)Block->ymin - EPSILON))
  2244.     {
  2245.       /* Test current block. */
  2246.  
  2247. #ifdef HFIELD_EXTRA_STATS
  2248.       Increase_Counter(stats[Ray_HField_Block_Tests_Succeeded]);
  2249. #endif
  2250.  
  2251.       if (dda_traversal(Ray, HField, nearP, &HField->Data->Block[z][x]))
  2252.       {
  2253.         if (HField->Type & IS_CHILD_OBJECT)
  2254.         {
  2255.           found = TRUE;
  2256.         }
  2257.         else
  2258.         {
  2259.           return(TRUE);
  2260.         }
  2261.       }
  2262.     }
  2263.  
  2264.     /* Step to next block. */
  2265.  
  2266.     x = nx;
  2267.     z = nz;
  2268.  
  2269.     Assign_Vector(nearP, farP);
  2270.  
  2271.     neary = fary;
  2272.   }
  2273.   while ((x >= 0) && (x < xmax) && (z >= 0) && (z < zmax) && (neary >= ymin) && (neary <= ymax));
  2274.  
  2275.   return(found);
  2276. }
  2277.  
  2278.